Accurate polarization within a unified Wannier function formalism 



in 
o 
o 

(N 
G 

m 

(N 



-a 
o 

(N 
> 
On 
OO 

m 

O 

in 
o 

-I— > 



i 

c 

o 
o 



X 



Massimiliano Stengel 1 and Nicola A. Spaldin 1 

1 Materials Department, University of California, Santa Barbara, CA 93106-5050, USA 

(Dated: February 2, 2008) 

We present an alternative formalism for calculating the maximally localized Wannier functions 
in crystalline solids, obtaining an expression which is extremely simple and general. In particular, 
our scheme is exactly invariant under Brillouin zone folding, and therefore it extends trivially to 
the F-point case. We study the convergence properties of the Wannier functions, their quadratic 
spread and centers as obtained by our simplified technique. We show how this convergence can be 
drastically improved by a simple and inexpensive "refinement" step, which allows for very efficient 
and accurate calculations of the polarization in zero external field. 

PACS numbers: 71.15.-m 



I. INTRODUCTION 



The representation of the one-particle electronic struc- 
ture of molecules and solids in terms of localized Wan- 
nieni orbitals is nowadays "enjoying a revival"— as a use- 
ful tool for many applicationo 3 ! 4 ! 5 ! 6 ! 7 ! 8 ! 9 . 10 ! 11 ! 12 ! 13 . The 
main impetus for this renewal of interest was given by 
the establishment, by King-Smith and Vanderbilt (KSV), 
of a formally exact relationship between the sum of the 
Wannier function (WF) centers and a gauge-invariant 
Berry phase, in the context of the modern theory of po- 
larization^. However, the intrinsic nonuniqueness in the 
Wannier function definition, and the difficulty in defining 
their centers within a periodic cell calculation, limited 
their practical use, until a particularly elegant method 
due to Marzari and Vanderbilt 3 (MV) became available 
some years ago. Their scheme allows one to obtain, in a 
given isolated or extended system, a unique set of maxi- 
mally localized Wannier functions that minimizes a well- 
defined spread functional. Moreover, the MV formal- 
ism provides as an important byproduct the positions of 
the WF centers, whose sum gives direct access to the 
macroscopic polarization of the physical system. The 
MV scheme, which became instantly popular, presents 
nevertheless an inconvenience, in that crystalline solids 
are treated on a different footing with respect to the case 
of, e.g., large disordered systems simulated at the T point 
only. The two prescriptions are indeed equivalent in the 
thermodynamic limit, but they formally differ when dis- 
crete Brillouin zone (BZ) samplings (or finite supercells 
for isolated objects) are used, which is necessarily the 
case in any practical calculation. In the first part of 
this work we show that there is nothing fundamental in 
this discrimination, i.e. that a given choice for a spread 
functional in T-sampled cells dictates unambiguously the 
mathematical expression in discrete fc-point space, and 
that invariance under "BZ folding" is the guideline which 
estabilishes the link. The resulting formalism is com- 
pletely general and, while being similar in spirit to the 
original (MV) one, presents a much simpler algebra. 

A more relevant issue affects more directly the mod- 
ern theory of polarization, and concerns the asymptotic 



convergence with respect to BZ sampling. It was already 
shown formally and numerically 3,15 that both methods 
(KSV and MV) for calculating the polarization of a 
molecule or a crystal in periodic boundary conditions 
(PBC) are plagued by a slow (D(L~ 2 ) convergence, where 
L is the linear dimension of the supercell containing the 
isolated molecule, or alternatively the resolution of the 
fc-point mesh. The problem is usually addressed, in the 
context of the Berry-phase KSV approach, by refining 
the fc-point grid along "stripes" in the Brillouin zone 
within a separate, non-selfconsistent calculation^. For 
finite systems, an extrapolation technique was recently 
proposed 16 , in which the 0(L~ 2 ) error is removed from 
the Wannier multipoles by performing a series of calcu- 
lations in cubic supercells of increasing size 3 ^. Both so- 
lutions are somewhat unsatisfactory, in that they require 
many calculations to be performed on the same system, 
with a cost that is higher than what is normally needed to 
converge total energies and densities. In the second part 
of this work we propose a simple "refining" procedure, 
which is able to provide an extremely accurate value for 
the center and spread once a well-localized set of maxi- 
mally localized Wannier functions is available. We show 
first formally and then by numerical examples that this 
technique, while requiring a very minor computational ef- 
fort, is able to outperform in terms of accuracy both the 
standard KSV Berry-phase approach and the alternative 
formula based on "unrefined" WF centers 3 ^. Finally, our 
derivation also provides a novel, intuitive interpretation 
of the position/localization operator in periodic bound- 
ary conditions and of its relationship to the correspond- 
ing well-known free-space operators. 



II. METHOD 

The theoretical basis for the MV approach rests on a 
continuum formulation, in which the space is infinitely 
extended in all directions; this translates to an infinitely 
dense Brillouin zone sampling in the case of crystalline 
solids. For practical calculations a finite sampling (or 
finite simulation supercell) is necessary, and MV give de- 
tailed prescriptions for the "discretization" of the rele- 
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vant mathematical objects (gradients and Laplacians in 
fc-space). We start here our alternative derivation from 
a slightly different viewpoint, i.e. we "discretize" the 
problem from the very beginning by choosing an appro- 
priate spread functional in the T-point case, and then 
work out the formulas in /c-space without making any 
further approximation. This approach leads automati- 
cally to a general and size-consistent formalism, that is 
invariant under BZ unfolding. 

We assume a Born- von Karman (BvK) supercell of vol- 
ume VbvK, which is a multiple of the primitive (P) unit 
cell of the crystal (of volume Vp) under study. In this 
system with periodic boundary conditions (PBC) there 
are N allowed Bloch vectors, where N is given by the 
ratio between the volumes: 



N 



BvK 
V P ' 



The generalized Bloch orbitals (which are not necessarily 
eigenstates of the Hamiltonian) are orthonormal on the 
primitive cell: 



C l k( r )V'nk(r)rfr = 5„ 



and can be written as usual: 



V'nk(r) 



Unk(r), 



where the u n w are periodic functions, and can be repre- 
sented on the reciprocal lattice of the P cell: 



Wnk(r) 



1 



E 



<G.r- 



u nk (G). 



G+k| 2 <B cut 



E C ut represents the plane- wave cutoff, while u n k(G) is 
the Fourier coefficient of the lattice-periodic part of the 
Bloch function: 



Mnk(G) = 



1 



-iG.: 



u nk (r)dr. 



We will use the BvK supercell for representing our Wan- 
nier functions: 



1 



%W = ^£?/W r ), 

k 

where the normalization constant is chosen so that these 
w n are orthonormal on the BvK supercell. A particularly 
simple relationship holds in reciprocal space: 



w n (G + k) = — !=u„ k (G). 

V N 



(1) 



We remind the reader that the reciprocal lattice of the 
BvK cell is spanned by all vectors of type b = G + k, 
which we will call b in the following, to distinguish them 
from the G vectors of the primitive reciprocal lattice. 

Eq. n does not define a unique set of Wannier func- 
tions, because of the gauge arbitrariness in the choice 



of the unitary representation of the Bloch vectors. This 
indeterminacy can be solved by defining a spread func- 
tional which depends explicitly on the gauge, so that 
the minimization of Q leads to a well defined set of local- 
ized orbitals with the desired properties. Berghold and 
coworkers** proposed a particularly simple and appeal- 
ing expression for Q and the related Wannier centers r„, 
which is valid for T-only BZ sampling in a lattice of gen- 
eral symmetry. Since our BvK supercell is sampled at 
r by construction, we can use the same expressions as 
Berghold, that in our notations read: 

r„ = ^Wibilmln^W (2a) 

i 

n = EE*- 2 ^-^!)- (2b) 



(i) 

Here z n are dimensionless complex numbers given by: 



z« = K|e ib *- r K 



and {bi,Wi} represents a small set of reciprocal lattice 
vectors b^ with weights w$. In the case of a cubic BvK 
supercell of edge L these quantities reduce to the i — 
1, ...3 primitive reciprocal-space vectors of the BvK cell 
and the weights are all equal: 



27T, 



(£)"■ 



while in the most general case of a triclinic cell a max- 
imum number of six independent vectors and weights 
must be used, according to the prescriptions given in 
Ref . and Hi 

With these notations and conventions in hand, we are 
now ready to write down a fc-space expression for r„ and 
VL. Both quantities depend directly on z„ , , and the key 
of the derivation is then the "Brillouin-zone unfolding" 
of this latter quantity. Using the same notation as MV: 

^mn bi) = (UmkKk+bJ, 

it is straightforward to derive a very simple expression 



for zl l) : 



= -Ym. 

N 4-' ' 



nn 



With this formula we can write our operational defini- 
tions of position and quadratic spread in fc-space: 



k,b 4 ) 
nn 



r n = -^w 4 bjmln [^J2 M « 

i k 



(3a) 
(3b) 



It is interesting to notice the strikingly close similarity 
between our expression for the centers (Eq. I3a|l and Eq. 
31 of MV, the only difference being the order in which 
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the complex logarithm and the average over fc-points is 
taker. 39 . We argue that the one proposed here is a more 
natural choice, since it retains the correct translational 
properties of their formula, while strictly enforcing size 
consistency. Size consistency means that the formalism 
gives mathematically identical answers for the fc-point 
representation and for the equivalent BvK real-space T- 
point representation. Our formula is correct by construc- 
tion, and extends exactly to the case of isolated systems 
with T-point sampling without any further algebra. 

Another advantage of our scheme is its simplicity, 
which becomes evident when taking the gradient of the 
spread functional with respect to an infinitesimal uni- 
tary rotation in a given k subspace. Thus, we consider 
thc transformation: 



(r) 



E 



'"'ink 



(rW (k) 

V / ran > 



where the rotation matrices U^ k ' are obtained by adding 
an infinitesimal antiHermitian matrix dW to the identity: 

C/( k > ~ l + dW (k) . 

The variation of the total spread with respect to this 
transformation is readily obtained in terms of the M^ k,b ^ 

(i) 

matrices and the phases of the z„ complex numbers: 

\dWWJmn N ^ l V mn °" + 



+Mt7 b " h ' ] *C^ ) - H.c 



(4) 



where H.c. stays for Hermitian-conjugate, and C„ is the 
phase: 



This expression for the gradient can easily be obtained 
by observing that, in Ea. I2bl one can write \z\ = ze~ 1 ^. 

We note that the for the spread functional (Eq. 12 b|) 
several possibilities exist, which are all equivalent in the 
thermodynamic limit^i. In the Appendix we briefly con- 
sider these alternatives, and we provide a formal deriva- 
tion of the gauge- invariant part of the spread^, which fur- 
ther evidences the close relationship of our formulation to 
the original MV scheme. Because of the exact mapping 
between the BvK supercell and the primitive one, we find 
it particularly natural to choose our Wannier functions 
to be real. Even if there is no formal proof that at the 
global minimum of f2 the Wannier functions are real, this 
is nevertheless a very reasonable assumption 3 , and allows 
one to fully take advantage of thc time-reversal symme- 
try, with significant gain in computational efficiency. 

(k) 

For the minimization of fl with respect to the Umn de- 
grees of freedom many efficient schemes are available 21 . 



We decided in this work to implement a damped dynam- 
ics algorithm, which allows for good control over the pro- 
cess, at the expense of requiring more human input for 
the optimal tuning of the two independent parameters 
(time step and friction). In antiferromagnetic MnO, a 
case that is known& to be difficult to converge, we were 
able to obtain this way a very accurate and symmetri- 
cal minimum (to machine precision) in a couple of thou- 
sand time steps, which required only a few minutes on a 
modern workstation. An even more appealing feature of 
the dynamical scheme is the availability of a mathemati- 
cally conserved constant of motion, which provides a very 
stringent test on the accuracy of the implementation. 



III. CONVERGENCE PROPERTIES 

Since the Wannier functions in an insulator are known 
to be exponentially localized in space 23 , similar conver- 
gence properties can be expected for any physical quan- 
tity that is extracted from this particular representation 
of the electronic structure. Instead, as we pointed out 
at the beginning, both the sum of Wannier centers and 
the Berry phase (which is formally related to the Wan- 
nier centers by the derivation in KSV) converge only as 
(D(L~ 2 ), and need special treatment whenever accurate 
values are needed. 

We will show in this section that this slow convergence 
is indeed not a intrinsic feature of the ground state elec- 
tronic structure of an extended system, and can be dra- 
matically improved by a simple, inexpensive and very 
general procedure. Before explaining our correction in 
detail, we will first provide an intuitive picture of the po- 
sition operator in PBC, which, as Resta showed 2 ^, is the 
"kernel" of both Berry phase and maximally localized 
Wannier function calculations. 

Let's consider a one-dimensional system of one single 
electronic state which we will assume to be well lo- 
calized within a periodic cell of length L. The expression 
for the fundamental, dimensionlcss complex number z is 
very similar to the 3D expression: 



(5) 



The average value of the position operator (Eq. I2a|) be- 
comes: 



L L 
x — — Immz = — d>, 
2tt 2tt 



(6) 



while the quadratic spread (Eq. I2b|) reduces to: 
L \2, 



By defining the charge density p{x) = \ip{x)\ 2 , it is easy 
to see that the following is trua^i: 



2tt 

p{x) sin[— [x — x)]dx = (7a) 
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FIG. 1: Pictorial representation of the position (red) and 
spread (black) operators as they are approximated by the 
Berry-phase formalism when working in periodic poundary 
conditions. 



n =(^r) 2 I p(x){2-2cos[^{x-x)]}dx (7b) 

These equations can be directly compared to the elemen- 
tary textbook definitions of the position and quadratic 
spread for a square-integrable electronic state in one di- 
mension (i.e. without PBC, the superscript F stays for 
"free-space" ) : 

/OO 
p(x)(x - x F )dx = (8a) 
-OO 

/OO 
p(x)(x-x F ) 2 dx (8b) 
-OO 

The resemblance is indeed striking, the only difference 
being the replacement of the x and x 2 operators with 
trigonometric functions that are periodic on the cell. 
This relationship between the (polynomial) free-space op- 
erators and the (trigonometric) PBC ones is made evi- 
dent in Fig. ^ where they are plotted together in order 
to show their close matching in a region surrounding the 
localized state. Indeed, by a Taylor expansion one ob- 
tains: 

/L\ 2 r „ „ ,2n ., 9 /27r\ 2 a: 4 
(-) [2-2cos( T .)]~^-2( T ) - + ... 

Thus, we arrived from a different starting point at the 
same (D(L~ 2 ) convergence of the position and spread, 
which has already been discussed in the literatur o 3 ! 15 . It 
is particularly clear from this derivation that the intrinsic 
property of the periodic lattice and the localized state are 



by no means responsible for the slow convergence, which 
is instead determined exclusively by the mathematical 
form of the PBC position operator. 

To end this section, we note that Eq. [7a] alone is not 
sufficient to define the center x, since also x + satisfies 
the same requirement. In the context of Equations I7al 
and I7bl a correct definition of x can be given as the 
points in the lattice which minimize Eq. I7bl (it is easy 
to show that this definition is identical to the standard 
one in Eq. [SJ) . Interestingly, from this point of view the 
position x can be thought of as an internal parameter of 
the formalism, which is implicit in the definition of the 
spread. 

IV. CORRECTION SCHEME 

Since we are working with Wannier functions which are 
expected to be well localized in space (as the ID state de- 
picted in Fig. there is actually no need to insist on 
using the PBC formulas for calculating Wannier centers. 
One could argue here that our "Wannier functions" are 
formally still periodic (although represented on a large 
BvK supercell), and since their Hilbert space is defined 
within PBC, only the action of PBC-allowed operators is 
justified on them. Actually, another point of view can be 
used. We recall that true Wannier functions are contin- 
uous functions in the full 3D (reciprocal) g-space. The 
mean value of a local operator V(r) in real space can be 
written 

(w n \V\w n ) = I \w n {Y)\ 2 V{r)dY (9) 

J All — space 

or, equivalently in g-space (p rl (q) is the continuous 
Fourier transform of the Wannier density p n ( r ) = 
K(r)| 2 ): 

(w n \V\w n ) = [ 7(q)*p n (q)dq (10) 

J\q\<Ec*t 

If the Wannier function is localized (i.e. zero beyond 
a given distance from its center) , the integral in Eq. [5] 
can be limited to a finite region of space, for example a 
cubic box centered around the region where the Wannier 
density is nonzero. The <?-space integral in Eq. I1UI can 
then be recast to a sum over a discrete set of reciprocal 
space vectors, which is also finite because of the plane- 
wave cut off, and the result is still exact. 

If the integration box is chosen to be smaller than the 
region where the Wannier density is nonzero, then the 
reciprocal-space sum carries an error which is due to the 
overlap between the tails of the Wannier functions and 
their (artificially) repeated images. This overlap depends 
on the decay properties of the localized state, and in par- 
ticular it goes exponentially to zero for increasing inte- 
gration box size whenever \w n ) is exponentially localized. 

In a standard DFT simulation of a periodic crystal, the 
discrete set of reciprocal-space Wannier function coeffi- 



5 



cients are denned by Eq.^ and converge to their thermo- 
dynamic limit as soon as the total charge density is con- 
verged. Then, the only effect of a further refinement of 
the fc-points mesh is an increase in the BvK cell volume, 
which leads to the progressive reduction of the overlap 
term discussed above. Therefore, assuming exponential 
decay for the Wannier functions, our technique allows for 
an exponential convergence of the calculated expectation 
value of any free-space operator. The natural "bounding 
box" for the integration domain in real space is, for a 
general lattice, a Wigner-Seitz BvK cell aligned on the 
Wannier center. With this choice, the discrete Fourier 
representation of a given local free-space operator (we 
use here again the standard conventions for normaliza- 
tions and Fourier transforms) is: 



V(h) 



1 



e- lhr V(r)dr, 



VBvK J Wigner-Seitz 

and the expectation value is simply given as 
K|V> n ) = Vb oJr £ V-*(b)p n (b) 

b 

Starting from a well-localized set of Wannier func- 
tions we can now define a "refined" spread operator 
fi' = tt' n , where the contribution from the individ- 
ual WF is: 

K= [ |r-r' n |V(r)dr. 

J Wigner—Seitz 

The 6-space expression for this formula can be derived 
starting from the Fourier series of a parabola in a one 
dimensional box: 



1 f~ .2nk 2 f L\ 

lJ_± cos{ — x)x dx = {^J 



L \ 2 2(-l) 



k 2 



(fc > 0) 



x dx = — , 



and is readily generalized to three dimensions using the 
same set of vectors and weights {bj,Wj} introduced in 
Section 2: 



ZVbvK WiRe 



i,fc>0 



2(-l) fe , 



febi 



k 2 



(11) 



The "refined" position f' n which appears in Eq.^Jis then 
again an internal parameter, which is defined by the min- 
imum of £l' n for a given p n (see the discussion at the end 
of Section 3) . By taking the gradient of d' n with respect 
to r' n one obtains that, at the minimum, the integral: 



Ar„ = 



(r-r' n )p n (r)dr 



vanishes. Consistently with the definition of the spread, 
this condition has to be enforced in reciprocal space, 
where this integral becomes: 



Ar„ = -2V Bv k ^2 Wjb, 



Re 



i,k>0 



.(-1)' 



e 4fcb '- F "p„(fcbO 



Wigner—Seitz 



(12) 

The stationary point can be obtained iteratively start- 
ing from a set of maximally localized Wannier functions 
and unrefined centers, by updating at every iteration f' n 
through the addition of Ar„ as calculated in Eg. 1 121 until 
convergence is reached. If the Wannier function is ex- 
actly zero in a region surrounding the boundary of the 
Wigner-Seitz cell, one iteration is sufficient to provide the 
exact value of the center, while for less converged cases 
up to ten iterations may be necessary to achieve machine 
precision. These iterations have anyway negligible cost, 
since the Fourier transform of the Wannier function on 
the BvK cell has to be evaluated only at the beginning 
(twice for each Wannier function to get the density in 
reciprocal space). 

Both expressions and f2' are in fact particular cases 
of a class of localization criteria which rely on individual 
Wannier densities only, through some generalized spread 
functional S: 

n 

A similar generalised, density-dependent spread can be 
used in practice to explore alternative localization crite- 
ria, like e.g. the maximal Coulomb self-repulsion of Ed- 
miston and Ruedenberg2£, or the orbital self-interaction 
as defined by Perdew and Zunger— . An article compar- 
ing such alternatives is under preparation^. 

Since the present free-space-like expressions for posi- 
tion and spread are more accurate than those derived in 
the first part of this work, one could wonder why we did 
not use them from the beginning. The reason is exclu- 
sively related to computational efficiency. In Eqns.|3]the 
localization algorithm involves operations on small J x J 
matrices only, where J is the number of bands in the 
primitive cell (the computationally intensive calculation 
of the Af( k,bi ) matrices has to be performed only once at 
the beginning of the iterative minimization). If, instead, 
the refined spread (or one of the alternative localization 
criteria discussed above) is used directly for localizing the 
Wannier functions, several Fourier transforms on the full 
Wannier (BvK) grid are required for each iteration, at 
a substantially higher cost. This expensive procedure is 
anyway not necessary for the scope of the present work, 
since the actual set of Wannier wavefunctions obtained 
from one localization method or the other coincide in 
practice to a high degree of accuracy 4 (in particular, the 
decay properties are expected to be very similar). There- 
fore, we find it most convenient to use this refinement 
step in a one-shot fashion once a set of maximally lo- 
calized Wannier functions is obtained within the more 
efficient localization functional (Eqns.[3J). 
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FIG. 2: Convergence of the Born effective charge Z* of oxygen 
in MgO as computed using the sum of unrefined Wannier 
centers (Eq. I3al circles) , the Berry-phase approach (Eq. 1131 
diamonds) and our new scheme (Eg. 1111 squares). The blow- 
up in the inset shows the improvement of our method with 
respect to the Berry-phase approach. 



FIG. 3: Convergence of the total spread in the case of MgO 
and Si. The improvement provided by the refined value 
(Eq. 1111 red line and squares) with respect to the standard 
"trigonometric" spread (Eq. I3bl black line and circles) is ap- 
parent. We provide for comparison the alternative spreads 
discussed in the Appendix, which are also based on the same 
trigonometric kernel and show similar, slow convergence prop- 
erties. 



NUMERICAL TESTS 



To demonstrate the effectiveness of our method we 
have chosen two examples which have been extensively 
studied in the literature: (i) the dynamical Born effec- 
tive charge of oxygen atoms in rocksalt MgO, and (ii) the 
spontaneous polarization of the ferroelectric, tetragonally 
distorted phase of KNDO3. Our calculations were per- 
formed within the local density approximation^, by us- 
ing norm-conserving Troullier and Martins^ pseudopo- 
tentials in the Kleinman and Bylander— form. A non- 
linear core correction^ was adopted for the Mg pseu- 
dopotential, while the K pseudopotential was generated 
in the 4s° ionized configuration with the semicore 3s, 3p 
orbitals included in the valence. We used the experi- 
mental lattice constants and atomic positions (ao=7.96 
a.u. for MgG^I, and the structural data for KNb03 from 
Ref. I32I) . We expanded the electronic ground state on a 
plane-wave basis up to a cutoff of 70 Ry. The BZ sam- 
pling was performed with T-centered simple cubic (or- 
thorombic) grids in reciprocal space for MgO (KNbOa), 
by taking into account the time-reversal symmetry only. 

We will compare the results as a function of fc-mesh 
resolution for three different methods for calculating the 
polarization: (i) the sum of Wannier centers as obtained 
by Eq. I3al (ii) the sum of refined Wannier centers as de- 
scribed in the previous section; (iii) the Berry-phase ap- 
proach. We note that the Berry-phase result can be read- 
ily obtained from the quantities that are already available 
in the localization formalism: 

r Berry = - ^ ^ w.b, ^ Im In det Af ( k ' b *) . (13) 

i k 



A. MgO 

The dynamical Born effective charge (Z*) of oxygen 
was calculated by the finite-difference method, i.e. by 
considering the difference in total polarization between 
the ideal centrosymmetric ground state and an atomic 
configuration where the oxygen sublattice was displaced 
by 1% of the cubic lattice constant along the x direc- 
tion. The atomic coordinates were prepared in such a 




0.36 b 1 1 1 , 1 1 1 , L 

3 4 5 6 7 

k mesh resolution 



FIG. 4: Convergence of the spontaneous polarization of the 
(tetragonal) ferroelectric phase of KNbOs, as computed us- 
ing the sum of unrefined Wannier centers (Eq. I3al circles), 
the Berry-phase approach (Eq. 1131 diamonds) and our new 
scheme (Eq. 1111 squares). 
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way that, in the ideal lattice, the O atom sits at the ori- 
gin, and in this case the electronic contribution to the 
polarization is exactly zero modulo a polarization quan- 
tum (all four Wannier functions are symmetrical about 
the O in this case). We compare in Fig. [2] the resulting 
value for Z* as calculated by the three different meth- 
ods (i-ih). The results clearly show that the sum of the 
unrefined Wannier centers can be very inaccurate, and 
even for the finest mesh the error is still large. MgO is 
probably a very unfortunate case in that each sp 3 -like 
Wannier function has a strongly asymmetric shape, and 
the errors in the individual centers do not cancel out ef- 
ficiently in the strained configuration, so that the total 
polarization carries an important deviation from the ex- 
act value. The Berry-phase calculation is a much better 
estimate, but in the inset it can be seen that the con- 
vergence is still relatively slow. As we explained in the 
preceding sections, it was already shown that the Berry- 
phase result converges only as 0(L~ 2 ), i.e. it shares the 
same asymptotic behaviour as the sum of the unrefined 
Wannier centers (albeit with a quite different prefactor 
in this particular case) . The sum of the refined Wannier 
centers instead shows an extremely fast convergence, and 
gives a very accurate result already for a 2 x 2 x 2 mesh. 
The value of Z* we obtain is -1.95, which is in excellent 
agreement with previous experimental and theoretical in- 
vestigations^. 

To complement our methodological test, we calculated 
also the refined value for the total quadratic spread as 
a function of /c-mesh resolution, and the results are re- 
ported in Fig. 21 It is clear that this quantity shares the 
same, excellent, convergence properties as the position 
operator (upper panel) . In the lower panel of the same 
figure we report for comparison the results of an analo- 
gous calculation of the total spread in bulk silicon. The 
convergence is slower than in MgO, as can be expected 
from the very different character of this covalent com- 
pound as compared to the highly ionic magnesium oxide, 
but the benefit that can be obtained through the use of 
the more accurate free-space definition of the spread is 
still very clear. The "unrefined" value of the spread is 
also compared to the alternative, very similar prescrip- 
tions discussed in the Appendix. 



B. KNb0 3 

We present in Fig.0|our results for the spontaneous po- 
larization of KNb03 . The sum of the unrefined Wannier 
centers is less inaccurate in this case, and is fairly close 
to the values obtained within the Berry phase formalism. 
The sum of the refined centers has, again, much better 
convergence properties than the two traditional methods. 
By increasing the mesh from 6x6x6 to 7x7x7 the 
value of the spontaneous polarization increases by 0.02 %, 
while within the same 7x7x7 mesh the traditional tech- 
niques carry an error which is two orders of magnitude 
higher. Extrapolating the 0(L~ 2 ) trend one can guess 



that ~ 70 fc-points along the reciprocal space stripes^* 
would be needed to achieve similar accuracy within the 
Berry-phase formalism. The final value we obtain, 0.38 
C/m 2 , compares very well with experimental data and 
previous theoretical investigation a 19 i 32 i 34 i 35 . 

VI. CONCLUSIONS 

In conclusion, we have derived a simple and general for- 
malism for the computation of maximally localized Wan- 
nier functions. We provide an intuitive picture of the con- 
vergence properties of this scheme and similar ones, relat- 
ing them to the Taylor expansion of elementary trigono- 
metric functions. We show that the convergence can be 
dramatically improved by a simple strategy based on the 
exponential localization of the Wannier functions in insu- 
lating materials. We expect our scheme to open the way 
to both accurate and efficient calculations of polarization 
properties in a wide range of physical systems, making 
the expensive linear-response approach or the relatively 
cumbersome non-self consistent calculation of "stripes" 
in reciprocal space unnecessary. 

As a final remark, we note that the Wannier function- 
based theory of polarization is becoming increasingly im- 
portant especially in disordered systems, where not only 
the global polarization but also the local bonding proper- 
ties and dipole moments may be interesting to follow dur- 
ing, e.g. a molecular dynamics simulation 5 . In these ap- 
plications the improved accuracy provided by our method 
could be an extremely valuable tool. 
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APPENDIX: DECOMPOSITION INTO 
INVARIANT, OFF-DIAGONAL AND DIAGONAL 
PARTS 

The form I2bl for the spread functional was chosen 
mainly because of its simplicity, and because it allows 
for a direct interpretation as the integral of the Wannier 
density multiplied by a real function on the BvK cell (see 
the discussion in Sec. 3). Unfortunately this expression 
does not lead to an elegant separation into invariant and 
non-invariant parts. However, this issue is readily solved 
by choosing an alternative definition of the spread: 

fW = (^) 2 (l-M 2 ), (A.l) 

which coincides with the T-point prescription of MV and 
which does allow for an exact separation of the invariant 
part. This choice still allows for the simple interpretation 



8 



based on cosine-like functions. If we define a function of 
x : 



2tt 

f(xo)= I p(x)cos[—(x - x )]dx 



(A.2) 



it is clear that when xq maximizes /, xq is automatically 
the Wannier center of Eq. [5a| Both expressions for the 
spread (SI as in Eq. |2H| and SImv discussed here) are 
consistent with the same value of xq at the minimum: 

fi= (^) 2 min2[l-/(x )] 



&MV — ( 7. I min[l - f 2 {x Q )} 

Moving on to 3D, the operational definition of the 
spread becomes: 



MV 



U) 

where it is easy to see that z n are nothing other than 
the matrix elements indicated as X nn , Y nn , Z nn in MV. 

Now, "folding" this expression in fc-space leads to a 
formula which is similar to Eq. 13 bl 



n 



Thinking in terms of the big BvK cell, this can be 
written equivalently as: 

= 4 Y,** NJ E l(R^|e- lbl ' r |R«)| 2 ), 



R., 



where the leading factor 1/N gives the spread per prim- 
itive cell, and the same notations as MV for the rt-th 
Wannier function at the R site, |Rn) are used. From 
this expression it is clear how to construct an obvious 
invariant quantity, f2/ (J is the number of bands in the 
primitive cell): 

Vi = ^J2^( NJ - E l(R«|e- lb " r |R'm)| 2 ), 



RR',fim 



and what remains to do is to "unfold" this formula in 
fc-space. A first simplification is trivial: 



i R,nm 



-ihi 



■|0m)| 2 ). 



A second simplification is obtained by reversing the for- 
mula between Eq. 5 and 6 of MV, leading to: 



R,nm k 

icitly | z | 2 = 2*2 and nol 



i R,nm k 

By writing explicitly \z\ 2 = 2*2 and noticing that 



we obtain the final expression in fc-space: 

^E^-E^EK'i 2 ) 

i mn k 



(A.3) 



which is exactly Eq. 34 of the MV paper. 

It is interesting to work out the remaining terms, fin 
and SIod, which are indicated as "diagonal" and "off- 
diagonal" parts in MV (we recall that VLd vanishes for 
a centrosymmetric system^). The easiest way is to first 
solve the expression for 



MV 



n 



D = X>( J -Ei< R ' 



n\e 



0n}| 2 ). 



R., 



By using an analogous algebra we readily arrive at the 
formula in fc-space: 



n 



MV 



^=E^-E^E|^ bi) | 2 ) 



from which it is very easy to evaluate Qqd'- 



OD 



e^e^eK^ 



N 

m^n k 



This means that Clj) is given by the following difference: 



1 



^=E^E^n' 



2 1 

~ N 



Eh 



(k,b s ) 



Comparing this formalism with the MV one, it is clear 
that 0/ and 91od are identical, while the terms f2 and 
Qd differ. This derivation provide in a certain sense a 
"unification" of the two, formerly distinct, MV prescrip- 
tions for the T-point case and in fc-space. The gradients 
with respect to unitary rotations of the Bloch orbitals 
are simply given by setting C — z (instead of C — z/\z\) 
in Eq. 4. 

To complete our discussion, we note that a third form 
of the one-particle quadratic spread was proposed by 
Resta and Sorella 3 ^, which leads to yet another opera- 
tional definition for the localization criterion: 



Mrs = ^ E E Wl ln 



The fc-space folding of this formula is straightforward, 
while the gradient is again given by Eq. 4, with C = 
z/\z\ 2 . All functionals fi, VLmv an d ^i?s are identical 
in the thermodynamic limit. For finite BvK cells they 
all share the same definition of the Wannier center. The 
resulting maximally localized Wannier functions them- 
selves are identical in cases where \z n \ are equal for all 
n = 1,..., J bands (e.g. bulk Si, centrosymmetric MgO 
crystal). The numerical value of the spread can differ 
slightly, because the higher orders in the Taylor expan- 
sion are different. Some examples concerning this dis- 
crepancy are reported in the main text (see, e.g., Fig.|3J). 
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